********************************************************************************
** Cederman, Galano, Girardin and Schvitz. War Did Make States.
** Article prepared for International Organization
** June 20, 2022
**
** Stata do-file: analysis_dyads.do
** Main script for dyadic analysis
** Required file paths set in runall.do
********************************************************************************

// Requires that dyad_prep_dyadic is run first such that dyad_cumul.dta and cumulneigh.dta intermediate files are ready to be merged in from the intermediate file directory


cd $ROOT
cd $INTERMEDIATEDIR
use "statedata.dta", clear 

xtset id year

bys year: egen n=count(area)
bys year: egen sumarea=sum(area)
bys year: egen marea=mean(area)

bys year: egen minarea=min(area)
bys year: egen maxarea=max(area)

bys year: egen rank = rank(-area)
gen pr = rank/(n+1)
gen lnpr = log(pr)


//////////////////////////////


// Compute frac measures and lognormal parameters
gen s = area/sumarea

bys year: egen sumpop = sum(population)

bys year: egen tconc = sum(s^2)
gen tfrac = 1 - tconc

bys year: egen sumpop2 = sum((population/sumpop)^2)
gen pfrac = 1 - sumpop2

bys year: egen sumlnarea = sum(lnarea)  
gen muhat = sumlnarea/n

gen lnrarea = lnarea/sumlnarea
bys year: egen lntconc = sum(lnrarea*lnrarea)


bys year: egen mlnarea = mean(lnarea)

bys year: egen sumsqdiff = sum((lnarea - muhat)^2) 
gen s2hat = sumsqdiff/n
gen shat = sqrt(s2hat)


sort year rank
list name area lnarea rank pr n muhat shat if year==1500 & rank<11, noobs
list name area lnarea rank pr n muhat shat if year==1790 & rank<11, noobs

xtset id year
sort id year

// define labels for regression output

label variable incinit1 "initiator"
label variable incattacked1 "target"
label variable inc1 "war"
label variable warXsize "war X state size"
label variable lngain "log gain"
label variable lnloss "log loss"
label variable lnloss_nd "log loss (nd)"
label variable llnumneighsa "num. neighbors"
// label variable llrthreata "neigh. threat "
label variable llarea "log state size"
label variable llage "log state age"
label variable llelevsd "log elev. SD"
label variable lcoastal "coastal access"
label variable llwargrowth "log cumul. war gains"
label variable llpeacegrowth "log cumul. peace gains"
label variable llwarshrink "log cumul. war losses"
label variable llpeaceshrink "log cumul. peace losses"
label variable warXwargrowth "war X cumul. war gains"
label variable warXpeacegrowth "war X cumul. peace gains"
label variable warXwarshrink "war X cumul. war losses"
label variable warXpeaceshrink "war X cumul. peace losses"
label variable peaceXwargrowth "peace X cumul. war gains"
label variable peaceXpeacegrowth "peace X cumul. peace gains"
label variable peaceXwarshrink "peace X cumul. war losses"
label variable peaceXpeaceshrink "peace X cumul. peace losses"
label variable lcentral "core"
label variable lurban "urbanization"


cd $ROOT
cd $OUTPUTDIR

	
	
// GENERATE TABLE 1: GROWTH AND LOSSES

eststo clear
global X1 "warXwargrowth peaceXwargrowth warXpeacegrowth peaceXpeacegrowth"
global X2 "warXwarshrink peaceXwarshrink warXpeaceshrink peaceXpeaceshrink"
global X1P "inc1##c.llwargrowth inc1##c.llpeacegrowth"
global X2P "inc1##c.llwarshrink inc1##c.llpeaceshrink"
global CONTROLS " llage lcoastal llnumneighsa llelevsd"
global CONTROLS2 "lcentral lurban"
global TIMEVARS1 "nogainyear ngspline*"
global TIMEVARS2 "nolossyear nlspline*"

reghdfe lngain  inc1 warXsize  llarea $CONTROLS $TIMEVARS1, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estimates store m1

reghdfe lngain  inc1 $X1 llarea $CONTROLS $TIMEVARS1, absorb(year) cluster(id)
reghdfe lngain  inc1 $X1 llarea $CONTROLS $TIMEVARS1, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estimates store m2

reghdfe lngain  inc1 $X1  llarea $CONTROLS $TIMEVARS1, absorb(year id) cluster(id)
reghdfe lngain  inc1 $X1  llarea $CONTROLS $TIMEVARS1, absorb(year id) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "Yes"
estimates store m3

reghdfe lnloss  inc1 warXsize  llarea $CONTROLS $TIMEVARS2, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estimates store m4

reghdfe lnloss  inc1 $X2 llarea $CONTROLS $TIMEVARS2, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estimates store m5

reghdfe lnloss  inc1 $X2 llarea $CONTROLS $TIMEVARS2, absorb(year id) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "Yes"
estimates store m6

cd $ROOT
cd $OUTPUTDIR
esttab m1 m2 m3 m4 m5 m6 using "statelevel.tex", replace /// 
b(4) se label pr2 star(+ 0.1 * 0.05 ** 0.01 *** 0.001) ///
nogaps eqlabels(none) obslast  depvars ///
title("Territorial change \label{tab:country}") ///
order (inc1 warXsize llarea $X1 $X2) /// $CONTROLS) ///
keep (inc1 warXsize llarea $X1 $X2) /// $CONTROLS) ///
 nonotes addnotes("Standard errors clustered on states in parentheses." /// 
 "$^{+}p<0.1$; $^{*}p<0.05$; $^{**}p<0.01.$; $^{***}p<0.001.$ ") ///
      scalars("geocontrols Geo. Controls" "yearfe Year FE" "statefe State FE")
	

// GENERATE FIGURE 9a,b
// INTERACTION PLOTS 

global PLOTOPTS "recast(line) recastci(rarea) title("") legend(order(1 "Peace" 2 "War")) plot1opts(lcolor(black) lpattern(dash)) plot2opts(lcolor(black))  ci1opts(color("102 205 0%50")) ci2opts(fcolor(red%80)) graphregion(color(white)) bgcolor(white)"

// GENERATE FIGURE 9a
reghdfe lngain  inc1##c.llarea $CONTROLS, absorb(year) cluster(id)
margins inc1, at (llarea = (3(1)15) llnumneighsa=2  llarea=6.3 llage=4.6 llelevsd=3.3 lcoastal=1)
marginsplot, xlab(3(1)15) $PLOTOPTS
cd $ROOT
cd $OUTPUTDIR
graph export country_area.png, replace

// GENERATE FIGURE 9b
reghdfe lngain  $X1P llarea $CONTROLS, absorb(year) cluster(id)
margins inc1, at (llwargrowth = (0(1)15) llnumneighsa=2 llarea=6.2 llage=4.6 llelevsd=3.3 lcoastal=1)
marginsplot, xlab(0(1)15) $PLOTOPTS
cd $ROOT
cd $OUTPUTDIR
graph export country_wargrowth.png, replace

reghdfe lnloss  inc1##c.llarea $CONTROLS, absorb(year) cluster(id)
margins inc1, at (llarea = (3(1)15) llnumneighsa=2  llage=4.6 llelevsd=3.3 lcoastal=1)
marginsplot, xlab(3(1)15) $PLOTOPTS
cd $ROOT
cd $OUTPUTDIR
graph export country_area_loss.png, replace

// GENERATE FIGURE 9d
reghdfe lnloss  llarea $X2P $CONTROLS, absorb(year) cluster(id)
margins inc1, at (llwarshrink = (0(1)13) llarea=6.3 llnumneighsa=2  llage=4.6 llelevsd=3.3 lcoastal=1)
marginsplot, xlab(0(1)13) $PLOTOPTS
cd $ROOT
cd $OUTPUTDIR
graph export country_llshrinkwar.png, replace 

////////////////////////////////////////////////////////////////////////////////////////////////////
// Partioning sample above and below median state size

// GEMERATE TABLE A14

eststo clear
global X1 "warXwargrowth peaceXwargrowth warXpeacegrowth peaceXpeacegrowth"
global CONTROLS " llage lcoastal llnumneighsa llelevsd"
global TIMEVARS1 "nogainyear ngspline*"


eststo clear


reghdfe lngain  inc1 warXsize  llarea $CONTROLS $TIMEVARS1 if abovemed==1, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estadd local part "Above Median"
estimates store m1

reghdfe lngain  inc1 warXsize  llarea $CONTROLS $TIMEVARS1 if abovemed==0, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estadd local part "Below Median"
estimates store m2

reghdfe lngain  inc1 $X1 llarea $CONTROLS $TIMEVARS1 if abovemed==1, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estadd local part "Above Median"
estimates store m3

reghdfe lngain  inc1 $X1 llarea $CONTROLS $TIMEVARS1 if abovemed==0, absorb(year) cluster(id)
estadd local geocontrols "Yes"
estadd local yearfe "Yes"
estadd local statefe "No"
estadd local part "Below Median"
estimates store m4


cd $ROOT
cd $OUTPUTDIR
esttab m1 m2 m3 m4 using "statelevel_median.tex", replace /// 
b(4) se label pr2 star(+ 0.1 * 0.05 ** 0.01 *** 0.001) ///
nogaps eqlabels(none) obslast  depvars ///
title("Territorial change with samples above and below median state size \label{tab:country:median}") ///
order (inc1 warXsize llarea $X1) /// $CONTROLS) ///
keep (inc1 warXsize llarea $X1) /// $CONTROLS) ///
 nonotes addnotes("Standard errors clustered on states in parentheses." /// 
 "$^{+}p<0.1$; $^{*}p<0.05$; $^{**}p<0.01.$; $^{***}p<0.001.$ ") ///
      scalars("geocontrols Geo. Controls" "yearfe Year FE" "statefe State FE" "part Subsample")


////////////////////////////////////////////////////////////////////////////////////////////////////
// Cox duration models
// 

gen larea = log(area)
gen lelevsd = log(elevationsd + 1)

capture drop warXsize 
gen warXsize = inc1*lnarea

gen lwarXsize = linc1*llarea
gen lwarXwargrowth = linc1*llwargrowth
gen lwarXpeacegrowth = linc1*llpeacegrowth
gen lwarXwarshrink = linc1*llwarshrink
gen lwarXpeaceshrink = linc1*llpeaceshrink




label variable linc1 "war"
label variable larea "log state size"
label variable warXsize "war X state size"
label variable lwargrowth "cumul. war gains"
label variable lpeacegrowth "cumul. peace gains"
label variable lwarshrink "log cumul. war losses"
label variable lpeaceshrink "log cumul. peace losses"

global CONTROLSD "lnumneighsa lcoastal lelevsd"

global XD "llwargrowth lwarXwargrowth llpeacegrowth lwarXpeacegrowth llwarshrink lwarXwarshrink llpeaceshrink lwarXpeaceshrink"


// GENERATE TABLE 2
eststo clear
stset year, failure(death)

stcox inc1 larea $CONTROLSD, nohr cluster(id) efron
estadd local geocontrols "Yes"
estimates store m1

stcox inc1 larea warXsize $CONTROLSD, nohr cluster(id) efron
estadd local geocontrols "Yes"
estimates store m2

stcox inc1 warXsize larea lwargrowth lpeacegrowth lwarshrink lpeaceshrink $CONTROLSD, nohr cluster(id) efron
estadd local geocontrols "Yes"
estimates store m3

esttab m1 m2 m3 using "cox_death.tex", replace  /// 
b(4) se label pr2 star(+ 0.1 * 0.05 ** 0.01 *** 0.001) ///
nogaps eqlabels(none) obslast nodepvars ///
title("Cox proportional hazard models of state death \label{tab:coxmodels}") ///
order (inc1 larea warXsize  lwargrowth lpeacegrowth lwarshrink lpeaceshrink) ///  $CONTROLSD) ///
keep (inc1 larea warXsize  lwargrowth lpeacegrowth lwarshrink lpeaceshrink) ///  $CONTROLSD) ///
 nonotes addnotes("Standard errors clustered on states in parentheses." /// 
 "$^{+}p<0.1$; $^{*}p<0.05$; $^{**}p<0.01.$; $^{***}p<0.001.$ ") ///
 scalars("geocontrols Geo. Controls")

 

// GENERATE FIGURE 10: INTERACTION PLOT
 
stcox inc1##c.larea $CONTROLSD, nohr strata(year)

margins inc1, at (larea = (3(1)15)  lnumneighsa=2 lelevsd=3.3) predict(xb)
marginsplot, xlab(3(1)15) $PLOTOPTS 
graph export cox_area_death.png, replace

stcox linc1##c.llarea linc1##c.llwargrowth linc1##c.llpeacegrowth linc1##c.llwarshrink linc1##c.llpeaceshrink  $CONTROLSD, nohr strata(year)

graph export cox_area_death.png, replace

********************************************************************************
********************************************************************************
** Robustness
********************************************************************************
********************************************************************************

** Add center/periphery and Urban population controls

// GENERATE TABLE A1

global CONTROLSD2 "lcentral lurban"


eststo clear
stset year, failure(death)

stcox inc1 larea $CONTROLSD $CONTROLSD2, nohr cluster(id) efron
estadd local geocontrols "Yes"
estimates store m1

stcox inc1 larea warXsize $CONTROLSD $CONTROLSD2, nohr cluster(id) efron
estadd local geocontrols "Yes"
estimates store m2

stcox inc1 warXsize larea lwargrowth lpeacegrowth lwarshrink lpeaceshrink $CONTROLSD $CONTROLSD2, nohr cluster(id) efron
estadd local geocontrols "Yes"
estimates store m3


esttab m1 m2 m3 using "cox_death_centper_urban_vars.tex", replace  /// 
b(4) se label pr2 star(+ 0.1 * 0.05 ** 0.01 *** 0.001) ///
nogaps eqlabels(none) obslast nodepvars ///
title("Cox proportional hazard models of state death \label{tab:coxmodels-cpu}") ///
order (inc1 larea warXsize lwargrowth lpeacegrowth lwarshrink lpeaceshrink lcentral lurban) ///  $CONTROLSD) ///
keep (inc1 larea warXsize  lwargrowth lpeacegrowth lwarshrink lpeaceshrink lcentral lurban) ///  $CONTROLSD) ///
 nonotes addnotes("Standard errors clustered on states in parentheses." /// 
 "$^{+}p<0.1$; $^{*}p<0.05$; $^{**}p<0.01.$; $^{***}p<0.001.$ ") ///
 scalars("geocontrols Geo. Controls")


********************************************************************************

// GENERATE FIGURE 2
// tfrac plot
tsline tconc if tconc>=0 & tconc<=0.3, yscale(range(0.0 0.3)) ylabel(0.0 (0.05) 0.3) graphregion(style(none) color(white)) ytitle("Territorial Concentration")
graph export herfindahl_plot_ab2.png, replace

// GENERATE FIGURE 3

twoway (kdensity lnarea if year==1590 & waryears<5, recast(area) range(0 20) lcolor(black) fcolor("102 205 0%5")) /// 
(kdensity lnarea if year==1690 & waryears<10, recast(area) range(0 20) lcolor(black) fcolor("102 205 0%20")) /// 
(kdensity lnarea if year==1790 & waryears<15, recast(area) range(0 20) lcolor(black) fcolor("102 205 0%45")) /// 
(kdensity lnarea if year==1590 & waryears>=5, recast(area) range(0 20) lcolor(black) fcolor("255 0 0%40")) ///
(kdensity lnarea if year==1690 & waryears>=10, recast(area) range(0 20) lcolor(black) fcolor("255 0 0%55")) ///
(kdensity lnarea if year==1790 & waryears>=15, recast(area) range(0 20) lcolor(black) fcolor("255 0 0%70")), ///
 legend( lab(1 "less war in 1590") lab(2 "less war in 1690") lab(3 "less war in 1790") lab(4 "more war in 1590") lab(5 "more war in 1690") lab(6 "more war in 1790") order(1 4 2 5 3 6)) ///
 graphregion(style(none) color(white)) ytitle(Density) xtitle(log state size)
graph export state_size_comparison_1500v1790_logscale.png, replace

